Setup

Libraries used

Load libraries necessary to run the script.

library(tidyverse) # general
## Warning: package 'dplyr' was built under R version 4.2.3
## Warning: package 'stringr' was built under R version 4.2.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggdist) # plot distributions
library(ggside) # plot distributions on the plot margins
## Registered S3 method overwritten by 'ggside':
##   method from   
##   +.gg   ggplot2
library(easystats) # estimate and process statistical models
## # Attaching packages: easystats 0.6.0 (red = needs update)
## ✖ bayestestR  0.13.1     ✔ correlation 0.8.4   
## ✖ datawizard  0.7.1      ✖ effectsize  0.8.3   
## ✖ insight     0.19.1.6   ✖ modelbased  0.8.6   
## ✖ performance 0.10.3     ✖ parameters  0.21.0  
## ✖ report      0.5.7      ✖ see         0.7.5   
## 
## Restart the R-Session and update packages in red with `easystats::easystats_update()`.
library(patchwork) # combine multiple plots
library(modelsummary) # for tables
## Warning: package 'modelsummary' was built under R version 4.2.3
## Warning in !is.null(rmarkdown::metadata$output) && rmarkdown::metadata$output
## %in% : 'length(x) = 2 > 1' in coercion to 'logical(1)'
## 
## Attaching package: 'modelsummary'
## 
## The following object is masked from 'package:parameters':
## 
##     supported_models
## 
## The following object is masked from 'package:insight':
## 
##     supported_models

options("modelsummary_format_numeric_latex" = "plain")
options(modelsummary_get = "broom")

library(see) # for plotting
library(ggraph) # needs to be loaded
library(ggdist) # plotting distributions
library(kableExtra) # export tables to latex
## 
## Attaching package: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library(correlation) # easy correlation

Load data



df <- read.csv("data/data_pruned_for_analysis_2024-04-26.csv")%>% 
  mutate(
    Typology = fct_relevel(video_primary_category, "urban", "green"),
    Valence = datawizard::rescale(Valence, to = c(1, 7), range = c(1, 9)),
    Arousal = datawizard::rescale(Arousal, to = c(1, 7), range = c(1, 9))
  ) %>% 
  rename("Participant" = Anonymised_ID ) 

length(unique(df$Participant))
## [1] 377

Analysis

video_summary <- df %>% 
   group_by(video_name_df) %>% 
  tally() %>% 
  arrange(n)  %>% 
  summarise(m = mean(n), s = sd(n), r = paste( range(n), collapse = " - "))

Each video was seen 23.16 times on average (SD = 1.876, range = 19 - 28).

Summary data per video (stimulus)

dfenv <- df %>% 
  group_by(video_name_df_keep, Typology) %>% 
  summarise_at(.vars = c( "Restoration", "WillingnessToWalk","Presence", "Width", "Structure", "Beauty","Interest", "Scenic", "Familiarity",  "Valence","Arousal", "Crowdedness" ), .funs = function(x) mean(x, na.rm = T)) %>% 
  ungroup() %>% 
  mutate(Typology = recode_factor(Typology, "urban" = "Urban", "green" = "Green" )) %>% 
  rename(`Willingness to walk` = WillingnessToWalk)

Participant

Here we only provide an overview to confirm we have loaded the correct amount of participants – as a sanity check. Participant characteristics are analysed and reported by a different file (Table 1).

dfsub <- df |>
  group_by(Participant) |>
  select(Participant, Age, Sex, Language, Crowd_Preference_Mean, starts_with("ipip"), SSA_Mean, SPS_Total_Mean, SIAS_Total_Mean, NSS, 
         Upbringing_Urban, Current_Urban) |>
  slice(1)

nrow(dfsub) # should be 377
## [1] 377

A total of 377 participants (Mean age = 30.8, SD = 10.0, range: [18, 67]; Sex: 49.3% females, 50.7% males, 0.0% other) participants were included for analysis.

(dfsub |>
  ggplot(aes(x = Age)) +
  geom_density(fill = "orange") +
  geom_vline(xintercept = mean(dfsub$Age), color = "red", linetype = "dashed") +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  theme_modern()) |
  (dfsub |>
    ggplot(aes(x = Sex, fill = Sex)) +
    geom_bar(stat = "count") +
    scale_y_continuous(expand = c(0, 0)) +
    scale_fill_see_d(guide = "none") +
    theme_modern())

Correlation plot



p_corrplot_u <- dfenv %>% 
  filter(Typology == "Urban") %>% 
  select(WLW = `Willingness to walk`, Crowdedness, Valence, Arousal, Restoration, Presence, Beauty, Interest, Structure, Width) %>%
  correlation(., bayesian = T) %>% 
  summary(redundant = FALSE) %>%
  plot() + scale_x_discrete(position = "top", expand = c(0, 0)) + labs(title = "Urban spaces") + scale_y_discrete(expand = c(0, 0)) +
  theme_modern() + theme(legend.position = "none") 
## Warning in genhypergeo_series_pos(U = c((n - 1)/2, (n - 1)/2), L = ((n + :
## Series not converged.

## Warning in genhypergeo_series_pos(U = c((n - 1)/2, (n - 1)/2), L = ((n + :
## Series not converged.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

p_corrplot_g <- dfenv %>% 
  filter(Typology == "Green") %>% 
  select(WLW = `Willingness to walk`, Crowdedness, Valence, Arousal, Restoration, Presence, Beauty, Interest, Structure, Width) %>%
  correlation(., bayesian = T) %>% 
  summary(redundant = FALSE) %>% 
  plot() + 
  scale_x_discrete(position = "top", expand = c(0, 0)) + labs(title = "Green spaces") + 
  scale_y_discrete(expand = c(0, 0)) +
  theme_modern() +
  theme(legend.position = "bottom") 
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

p_corr <- p_corrplot_u / p_corrplot_g + plot_annotation(tag_levels = "A", title = "Correlation Matrix of emotions and environmental appraisals", subtitle = "Self-reported measures", caption = "Note: WLW = willingness to walk.")
p_corr


ggsave(p_corr, filename = "output/plots_corr.jpg", width = 13, height = 10, scale = 1, dpi = 300)

The environments

Helper function

Calculates bayesian t-test and reports BF, to be used when making tables of descriptives.


custom_ttest <- function(x, ...) {
    # Construct vectors of data y, and groups (strata) g
    y <- unlist(x)
    g <- factor(rep(1:length(x), times=sapply(x, length)))
    if (is.numeric(y)) {
        # For numeric variables, perform a standard 2-sample t-test
        # p <- t.test(y ~ g)$p.value
        p <- BayesFactor::ttestBF(x=y, y= g)
        p <- p@bayesFactor$bf
            # print(p)
    } else {
        # For categorical variables, perform a chi-squared test of independence
        p <- chisq.test(table(y, g))$p.value
    }
    # Format the p-value, using an HTML entity for the less-than sign.
    # The initial empty string places the output on the line below the variable label.
    c("", sub("<", "&lt;", format.pval(p, digits=3, eps=0.001)))
}

count_missing <- function(x) {
  sum(is.na(x))
}

order_names <- c( "WillingnessToWalk",
                  "Arousal",
                  "Valence",
                  "Restoration", 
                  "Presence" ,
                  "Beauty",
                  "Crowdedness",
                  "Width",
                  "Structure",
                  "Scenic",
                  "Interest",
                  "Familiarity")

ttest_bf <- df %>% 
  select(video_name, Typology, Restoration:Arousal) %>% 
  pivot_longer(cols = 3:14, values_transform = as.numeric) %>% 
  filter(!is.na(value)) %>% 
  group_by(Typology, name) %>% 
  nest() %>% 
  spread(key = Typology, value = data) %>% 
  mutate(
    t_test = map2(urban, green, ~{BayesFactor::ttestBF(.x$value, .y$value)@bayesFactor}[,-3]),
    urban = map(urban, nrow),
    green = map(green, nrow)
  ) %>% 
  unnest(cols = c(urban, green, t_test)) %>% 
  arrange(factor(name, levels =order_names)) %>% 
  ungroup() %>% 
  select(`BF*` = bf)
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.


length(ttest_bf)
## [1] 1
length(3:14)
## [1] 12

stimuli = df %>% group_by(video_name) %>% slice(1)%>% ungroup() %>% group_by(Typology) %>%count() 

df %>% 
  select( Typology, order_names) %>% 
  mutate(Typology = recode(Typology, "urban" = "Urban (397 stimuli)", "green"= "Green (83 stimuli)")) %>% 
  as.data.frame() %>% 
  datasummary(formula = All(.)  ~ Typology * ( (`Unique (\\#)` = N) +  Mean + SD ), # (`Missing (%)` = PercentMissing) +
              add_columns = ttest_bf, 
              title = "Descriptive statistics of emotional reactions and environmental appraisals of the stimuli, grouped by typology. \\label{tab:responseDescriptives}",
              notes = '(*) Bayes Factor of a bayesian t-test',
              output = 'output/table_response_descriptives.tex')
## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
##   # Was:
##   data %>% select(order_names)
## 
##   # Now:
##   data %>% select(all_of(order_names))
## 
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

plots


penv <- dfenv %>% 
  pivot_longer(Restoration : Crowdedness) %>% 
  ggplot(aes(x = name, y = value, colour = Typology)) +
  # geom_boxplot( position = position_dodge()) +
  # geom_line(aes(group = name_df_keep),  alpha =.2)+
  geom_jitter(position = position_jitterdodge(), alpha = 0.5) +
  ggdist::stat_halfeye(position = position_dodge(), aes(fill = Typology),
    ## custom bandwidth
    adjust = .5,
    ## adjust height
    width = .6,
    ## move geom to the right
    justification = -.4,
    ## remove slab interval
    .width = 0,
  point_colour = NA) +
  scale_color_manual(values = c("Urban" = col_urban, "Green" = col_green)) +
  scale_fill_manual(values = c("Urban" = col_urban, "Green" = col_green)) +
  coord_flip() +
  scale_y_continuous(breaks = 1:7) +
  geom_hline(yintercept = 4, linetype = "dashed") +
  # scale_y_discrete(position = "top") + 
  theme_modern() +
  theme(legend.position = "bottom") +
  labs(title="Average rating per stimulus", y = "Rating", x = "Dimension")

penv

Joint plot

pp <- (p_corrplot_u / p_corrplot_g)

pp_penv <- penv + pp
pp_penv[[2]] <- pp_penv[[2]] + plot_layout(tag_level = 'new')

pp_penv <- pp_penv + plot_annotation(tag_levels = c('A', '1'), 
                            title = "Environmental appraisals and emotions", 
                            subtitle = "Self-reported measures", caption = "Note: WLW = willingness to walk.")

ggsave(pp_penv, filename = "output/plots_envstats.jpg", width = 13, height = 13/2, dpi = 300, scale = 1.5)

Determinants of Subjective Crowdedness

Models

try(rm(preds, params, mods, vars), silent = T) # remove objects from previous runs, to be safe.
## Warning in rm(preds, params, mods, vars): object 'preds' not found
## Warning in rm(preds, params, mods, vars): object 'params' not found
## Warning in rm(preds, params, mods, vars): object 'mods' not found
## Warning in rm(preds, params, mods, vars): object 'vars' not found
preds <- data.frame()
params <- data.frame()
mods <- list()
outcome_var <- "Crowdedness"
vars <- c("Beauty", "Scenic","Width", "Structure") # replaced familiarity with Beauty
i=1
for (i in 1:length(vars)) {
  var <- vars[i]
  print(var)

  f <- paste0(outcome_var, " ~ log1p(video_mean_pedcounts) * ", var, " +  log1p(video_mean_pedcounts) * Typology + ", var, "* Typology + Familiarity  + (", var, "|Participant)")

  model <- lme4::lmer(as.formula(f), data = df )
  # model4 <- ordinal::clmm(as.formula( f ), data = df %>% mutate(Crowdedness = as.factor(Crowdedness)))
  mods[[i]] <- model
  
  # Parameters
  param <- parameters::parameters(model, effects = "fixed")
  param$Variable <- var
  params <- rbind(params, param)

  # Predicted
  pred <- estimate_relation(model, at = c("video_mean_pedcounts", var, "Typology"), length = c(50, 7))
  names(pred)[names(pred) == var] <- "Value"
  pred$Variable <- var
  pred$Label <- paste0("Main effect: ", format_p(param$p[3]), "\nInteraction: ", format_p(param$p[4]))

  preds <- rbind(preds, pred)
}
## [1] "Beauty"
## [1] "Scenic"
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00436804 (tol = 0.002, component 1)
## [1] "Width"
## [1] "Structure"
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00249185 (tol = 0.002, component 1)

preds_crowd <- preds

plots

rm(p_crowd)
## Warning in rm(p_crowd): object 'p_crowd' not found
p_crowd <- preds_crowd |>
  # filter(Predicted < 7) |>
  mutate(Value = as.factor(Value)) |>
  ggplot(aes(x = video_mean_pedcounts, y = Predicted)) +
  stat_density_2d(data = df, aes(y = Crowdedness, fill = ..density..), geom = "raster", contour = FALSE) +
  # geom_ribbon(aes(ymin = CI_low, ymax = CI_high, group = interaction(Typology, Value)), alpha =.2) +
  geom_line(aes(color = Typology, size = Value, group = interaction(Typology, Value))) +
  scale_fill_gradientn(colours = c("white", "grey"), guide = "none") +
  scale_size_manual(values = seq(from = .05, to = 2, length.out= 7), breaks= c(1:7)) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0), breaks = c(0, 1, 3, 5, 7)) +
  labs(x = "Number of Pedestrians", y = "Perceived Crowdedness") +
  labs(alpha = "Rating", title = "Predictors of Crowdedness") +
  theme_modern(axis.title.space = 10) +
  theme(
    strip.placement = "outside",
    strip.text.x = element_text(hjust = 0, size = 12),
    axis.title.y = element_text(face = "bold"),
    plot.title = element_text(face = "bold", hjust = 0.5)
  ) +
  ggnewscale::new_scale_fill() +
  scale_color_manual(values = c("urban" = col_urban, "green" = col_green )) + # col_green
  scale_fill_manual(values = c("urban" = col_urban, "green" = col_green)) +
  ggside::geom_xsidedensity(data = df, aes(y = after_stat(scaled), fill = Typology), color = NA, alpha = 0.5, adjust = 2) +
  ggside::geom_ysidedensity(data = df |>
    select(all_of(c("Typology", vars))) |>
    pivot_longer(-Typology, values_to = "Predicted", names_to = "Variable"), aes(x = after_stat(scaled), fill = Typology), color = NA, alpha = 0.5, adjust = 2) +
  ggside::theme_ggside_void() +
  facet_wrap(~Variable, strip.position = "top", ncol = 2)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

p_crowd
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Removed 12 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 22 rows containing non-finite values (`stat_density()`).

tables


names(mods) <- vars
names(mods) <- paste("Model",1:length(mods))

pretty_vars <- map(mods, .f = format_parameters) %>% unname() %>% unlist() 
pretty_vars
##                                       (Intercept) 
##                                     "(Intercept)" 
##                       log1p(video_mean_pedcounts) 
##                    "video mean pedcounts [log1p]" 
##                                            Beauty 
##                                          "Beauty" 
##                                     Typologygreen 
##                                "Typology [green]" 
##                                       Familiarity 
##                                     "Familiarity" 
##                log1p(video_mean_pedcounts):Beauty 
##           "video mean pedcounts [log1p] × Beauty" 
##         log1p(video_mean_pedcounts):Typologygreen 
## "video mean pedcounts [log1p] × Typology [green]" 
##                              Beauty:Typologygreen 
##                       "Beauty × Typology [green]" 
##                                       (Intercept) 
##                                     "(Intercept)" 
##                       log1p(video_mean_pedcounts) 
##                    "video mean pedcounts [log1p]" 
##                                            Scenic 
##                                          "Scenic" 
##                                     Typologygreen 
##                                "Typology [green]" 
##                                       Familiarity 
##                                     "Familiarity" 
##                log1p(video_mean_pedcounts):Scenic 
##           "video mean pedcounts [log1p] × Scenic" 
##         log1p(video_mean_pedcounts):Typologygreen 
## "video mean pedcounts [log1p] × Typology [green]" 
##                              Scenic:Typologygreen 
##                       "Scenic × Typology [green]" 
##                                       (Intercept) 
##                                     "(Intercept)" 
##                       log1p(video_mean_pedcounts) 
##                    "video mean pedcounts [log1p]" 
##                                             Width 
##                                           "Width" 
##                                     Typologygreen 
##                                "Typology [green]" 
##                                       Familiarity 
##                                     "Familiarity" 
##                 log1p(video_mean_pedcounts):Width 
##            "video mean pedcounts [log1p] × Width" 
##         log1p(video_mean_pedcounts):Typologygreen 
## "video mean pedcounts [log1p] × Typology [green]" 
##                               Width:Typologygreen 
##                        "Width × Typology [green]" 
##                                       (Intercept) 
##                                     "(Intercept)" 
##                       log1p(video_mean_pedcounts) 
##                    "video mean pedcounts [log1p]" 
##                                         Structure 
##                                       "Structure" 
##                                     Typologygreen 
##                                "Typology [green]" 
##                                       Familiarity 
##                                     "Familiarity" 
##             log1p(video_mean_pedcounts):Structure 
##        "video mean pedcounts [log1p] × Structure" 
##         log1p(video_mean_pedcounts):Typologygreen 
## "video mean pedcounts [log1p] × Typology [green]" 
##                           Structure:Typologygreen 
##                    "Structure × Typology [green]"


md <- modelsummary(mods, coef_omit = "Cor|SD", coef_rename = rename_explanator,
             statistic = "[{conf.low}, {conf.high}]",
             estimate = "{estimate}{stars} ",
             shape = term ~ model + statistic,
             gof_map = c("nobs", "r2.conditional","r2.marginal"), fmt = 3, 
             output = "dataframe",
             ) 
md  
md <- md %>% 
  select(-part) %>% 
  mutate(term = if_else(term %in% c("Num.Obs.", "R2 Cond.", "R2 Marg."), term, pretty_vars[md$term] %>% unname())) %>% 
  mutate(term = str_replace_all(term, "video mean pedcounts", "Ped.Counts")) %>% 
  mutate(term = str_replace_all(term, "Crowdedness", "Subj. Crowding")) 
md

md %>% 
  kable(format = "latex", caption = 'Linear mixed-effects models of Subjective Crowding.', label = "models_perceived_crowds",
        booktabs = T, col.names = c("Term", rep(c("Est.", "95% CI"), 4))) %>% 
  add_header_above(setNames(c("", 2, 2, 2, 2), c("", names(mods))) ) %>% 
  kable_classic_2(latex_options = c("striped", "scale_down")) %>%
  kableExtra::row_spec(., row = nrow(md)-3, hline_after = TRUE) %>% 
  # landscape() %>% 
  save_kable(file = "output/table_models_perceived_crowds.tex")


# to plot in the html version
modelsummary(mods, stars = TRUE) 
Model 1  Model 2  Model 3  Model 4
(Intercept) 1.534*** 1.440*** 1.791*** 3.706***
(0.093) (0.088) (0.087) (0.102)
log1p(video_mean_pedcounts) 1.889*** 1.920*** 1.864*** 1.199***
(0.046) (0.043) (0.040) (0.045)
Beauty 0.007
(0.020)
Typologygreen −0.174 −0.112 −0.148 −0.798***
(0.144) (0.136) (0.102) (0.130)
Familiarity 0.029*** 0.026*** 0.041*** 0.057***
(0.007) (0.007) (0.007) (0.007)
log1p(video_mean_pedcounts) × Beauty −0.045***
(0.010)
log1p(video_mean_pedcounts) × Typologygreen 0.616*** 0.639*** 0.629*** 0.651***
(0.052) (0.052) (0.049) (0.048)
Beauty × Typologygreen −0.063**
(0.024)
Scenic 0.035+
(0.019)
log1p(video_mean_pedcounts) × Scenic −0.056***
(0.010)
Scenic × Typologygreen −0.083***
(0.023)
Width −0.052**
(0.018)
log1p(video_mean_pedcounts) × Width −0.052***
(0.009)
Width × Typologygreen −0.074***
(0.018)
Structure −0.398***
(0.020)
log1p(video_mean_pedcounts) × Structure 0.040***
(0.009)
Structure × Typologygreen 0.060**
(0.021)
SD (Intercept Participant) 0.767 0.762 0.684 0.692
SD (Beauty Participant) 0.152
Cor (Intercept~Beauty Participant) −0.844
SD (Scenic Participant) 0.156
Cor (Intercept~Scenic Participant) −0.848
SD (Width Participant) 0.143
Cor (Intercept~Width Participant) −0.786
SD (Structure Participant) 0.156
Cor (Intercept~Structure Participant) −0.828
SD (Observations) 1.132 1.129 1.106 1.038
Num.Obs. 11107 11104 11105 11102
R2 Marg. 0.565 0.564 0.575 0.620
R2 Cond. 0.631 0.632 0.647 0.688
AIC 35221.3 35171.4 34766.8 33384.3
BIC 35309.1 35259.2 34854.6 33472.1
ICC 0.2 0.2 0.2 0.2
RMSE 1.11 1.10 1.08 1.01
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Determinants of Beauty

rm(preds, mods, vars)
preds <- data.frame()
dat <- data.frame()
mods <- list()
vars <- c("Crowdedness", "Structure", "Scenic", "Width") #  "Interest",
i = 1
for (i in 1:length(vars)) {
  var <- vars[i]
  print(var)

  f <- paste0("Beauty ~ poly(", var, ", 2) * Typology + Familiarity + (", var, "|Participant)")

  model <- lme4::lmer(as.formula(f), data = na.omit(df[c(var, "Beauty", "Typology", "Familiarity", "Participant")]))
  summary(model)
  
  mods[[i]] <- model

  # Parameters
  param <- parameters::parameters(model, effects = "fixed")

  # Predicted
  pred <- estimate_relation(model, at = c(var, "Typology"), length = c(30))
  names(pred)[names(pred) == var] <- "Value"
  pred$Variable <- var
  pred$Label <- paste0(
    "Objective Crowding: ", format_p(param$p[2]),
    "\n", var, ": ", format_p(param$p[3]),
    "\nInteraction: ", format_p(param$p[4])
  )

  # Data for density
  dat <- rbind(dat, data.frame(Value = df[[var]], Predicted = df$Beauty, Variable = var))
  preds <- rbind(preds, pred)
}
## [1] "Crowdedness"
## [1] "Structure"
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0043363 (tol = 0.002, component 1)
## [1] "Scenic"
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0198563 (tol = 0.002, component 1)
## [1] "Width"
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0109548 (tol = 0.002, component 1)

preds_beauty <- preds

plots


p_beauty <- preds_beauty |>
  ggplot(aes(x = Value, y = Predicted)) +
  stat_density_2d(data = dat, aes(fill = ..density..), geom = "raster", contour = FALSE) +
  geom_line(aes(color = Typology), size = 1) +
  scale_fill_gradientn(colours = c("white", "grey"), guide = "none") +
  scale_color_manual(values = c("urban" = col_urban, "green" = col_green)) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  labs(y = "Beauty", title = "Predictors of Beauty") +
  theme_modern(axis.title.space = 10) +
  theme(
    axis.title.x = element_blank(),
    strip.placement = "outside",
    strip.text.x = element_text(hjust = 0.5, size = 12, face = "plain"),
    axis.title.y = element_text(face = "bold"),
    plot.title = element_text(face = "bold", hjust = 0.5)
  ) +
  ggnewscale::new_scale_fill() +
  scale_fill_manual(values = c("urban" = col_urban, "green" = col_green)) +
  ggside::geom_ysidedensity(data = mutate(df, Predicted = Beauty), aes(x = after_stat(scaled), fill = Typology), color = NA, alpha = 0.5, adjust = 2) +
  ggside::geom_xsidedensity(data = df |>
    select(all_of(c("Typology", vars))) |>
    pivot_longer(-Typology, values_to = "Value", names_to = "Variable"), aes(y = after_stat(scaled), fill = Typology), color = NA, alpha = 0.5, adjust = 2) +
  ggside::theme_ggside_void() +
  facet_wrap(~Variable, strip.position = "top", scales = "free_x")

p_beauty
## Warning: Removed 34 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 12 rows containing non-finite values (`stat_density()`).
## Warning: Removed 22 rows containing non-finite values (`stat_density()`).

Tables

names(mods) <- vars
names(mods) <- paste("Model", 1:length(mods))

pretty_vars <- map(mods, .f = format_parameters) %>% unname() %>% unlist() 
pretty_vars
##                                   (Intercept) 
##                                 "(Intercept)" 
##                         poly(Crowdedness, 2)1 
##                    "Crowdedness [1st degree]" 
##                         poly(Crowdedness, 2)2 
##                    "Crowdedness [2nd degree]" 
##                                 Typologygreen 
##                            "Typology [green]" 
##                                   Familiarity 
##                                 "Familiarity" 
##           poly(Crowdedness, 2)1:Typologygreen 
## "Crowdedness [1st degree] × Typology [green]" 
##           poly(Crowdedness, 2)2:Typologygreen 
## "Crowdedness [2nd degree] × Typology [green]" 
##                                   (Intercept) 
##                                 "(Intercept)" 
##                           poly(Structure, 2)1 
##                      "Structure [1st degree]" 
##                           poly(Structure, 2)2 
##                      "Structure [2nd degree]" 
##                                 Typologygreen 
##                            "Typology [green]" 
##                                   Familiarity 
##                                 "Familiarity" 
##             poly(Structure, 2)1:Typologygreen 
##   "Structure [1st degree] × Typology [green]" 
##             poly(Structure, 2)2:Typologygreen 
##   "Structure [2nd degree] × Typology [green]" 
##                                   (Intercept) 
##                                 "(Intercept)" 
##                              poly(Scenic, 2)1 
##                         "Scenic [1st degree]" 
##                              poly(Scenic, 2)2 
##                         "Scenic [2nd degree]" 
##                                 Typologygreen 
##                            "Typology [green]" 
##                                   Familiarity 
##                                 "Familiarity" 
##                poly(Scenic, 2)1:Typologygreen 
##      "Scenic [1st degree] × Typology [green]" 
##                poly(Scenic, 2)2:Typologygreen 
##      "Scenic [2nd degree] × Typology [green]" 
##                                   (Intercept) 
##                                 "(Intercept)" 
##                               poly(Width, 2)1 
##                          "Width [1st degree]" 
##                               poly(Width, 2)2 
##                          "Width [2nd degree]" 
##                                 Typologygreen 
##                            "Typology [green]" 
##                                   Familiarity 
##                                 "Familiarity" 
##                 poly(Width, 2)1:Typologygreen 
##       "Width [1st degree] × Typology [green]" 
##                 poly(Width, 2)2:Typologygreen 
##       "Width [2nd degree] × Typology [green]"

md <- modelsummary(mods, coef_omit = "Cor|SD", coef_rename = rename_explanator,
             statistic = "[{conf.low}, {conf.high}]",
             estimate = "{estimate}{stars} ",
             shape = term ~ model + statistic,
             gof_map = c("nobs", "r2.conditional","r2.marginal"), fmt = 3, 
             output = "dataframe",
             ) 
md  
md <- md %>% 
  select(-part) %>% 
  mutate(term = if_else(term %in% c("Num.Obs.", "R2 Cond.", "R2 Marg."), term, pretty_vars[md$term] %>% unname())) %>% 
  mutate(term = str_replace_all(term, "video mean pedcounts", "Ped.Counts")) %>% 
  mutate(term = str_replace_all(term, "Crowdedness", "Subj. Crowding")) 


md
md %>%   
   kable(format = "latex", 
        caption = 'Non-linear mixed effects models of Perceived Beauty.', 
        label = "models_beauty",
        booktabs = T, col.names = c("Term", rep(c("Est.", "95% CI"), 4))) %>% 
  add_header_above(setNames(c(0, 2, 2, 2, 2), c("", paste("Model", 1:length(mods)))) ) %>% 
  kable_classic_2(latex_options = c("striped", "scale_down")) %>%
  kableExtra::row_spec(., row = nrow(md)-3, hline_after = TRUE) %>% 
  # landscape() %>% 
  save_kable(file = "output/table_models_perceived_beauty.tex")



modelsummary(mods, stars = TRUE)
Model 1  Model 2  Model 3  Model 4
(Intercept) 3.828*** 3.981*** 4.462*** 4.047***
(0.049) (0.045) (0.024) (0.045)
poly(Crowdedness, 2)1 5.688*
(2.219)
poly(Crowdedness, 2)2 −21.799***
(1.593)
Typologygreen 1.321*** 1.044*** 0.240*** 1.154***
(0.040) (0.038) (0.025) (0.034)
Familiarity 0.130*** 0.114*** 0.033*** 0.100***
(0.008) (0.008) (0.005) (0.008)
poly(Crowdedness, 2)1 × Typologygreen −25.855***
(4.278)
poly(Crowdedness, 2)2 × Typologygreen 11.589**
(3.870)
poly(Structure, 2)1 53.159***
(2.181)
poly(Structure, 2)2 3.441*
(1.564)
poly(Structure, 2)1 × Typologygreen 3.927
(4.626)
poly(Structure, 2)2 × Typologygreen −5.585
(4.005)
poly(Scenic, 2)1 136.398***
(1.275)
poly(Scenic, 2)2 −2.354*
(0.966)
poly(Scenic, 2)1 × Typologygreen 3.271
(2.974)
poly(Scenic, 2)2 × Typologygreen −4.676+
(2.469)
poly(Width, 2)1 58.481***
(2.308)
poly(Width, 2)2 4.007**
(1.555)
poly(Width, 2)1 × Typologygreen −17.420***
(3.424)
poly(Width, 2)2 × Typologygreen 2.241
(3.319)
SD (Intercept Participant) 0.756 0.997 0.536 1.042
SD (Crowdedness Participant) 0.149
Cor (Intercept~Crowdedness Participant) −0.584
SD (Structure Participant) 0.167
Cor (Intercept~Structure Participant) −0.845
SD (Scenic Participant) 0.091
Cor (Intercept~Scenic Participant) −0.915
SD (Width Participant) 0.177
Cor (Intercept~Width Participant) −0.856
SD (Observations) 1.279 1.215 0.727 1.198
Num.Obs. 11107 11102 11104 11105
R2 Marg. 0.151 0.195 0.688 0.195
R2 Cond. 0.371 0.519 0.798 0.542
AIC 38092.8 36871.8 25237.6 36647.2
BIC 38173.2 36952.2 25318.1 36727.7
ICC 0.3 0.4 0.4 0.4
RMSE 1.25 1.19 0.71 1.17
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Determinants of Psychological Reaction

rm(preds, dat, mods, vars)
preds <- data.frame()
dat <- data.frame()

mods <- list()
outcome_vars <- c("WillingnessToWalk", "Valence", "Arousal", "Restoration",  "Presence")
i=1
for (i in 1:length(outcome_vars)) {
  resp <- outcome_vars[i]
  print(resp)

  f <- paste(resp, " ~ Typology / (poly(Crowdedness, 2) * Beauty) + Beauty + Familiarity + (1|Participant)")
  f
  model <- lmerTest::lmer(as.formula(f), data = na.omit(df[c(resp, "Beauty", "Typology", "Crowdedness", "Familiarity", "Participant")]))
  mods[[i]] <- model
  param <- parameters::parameters(model, effects = "fixed")
  param
  # Predictions
  preds <- estimate_relation(model, at = c("Typology", "Crowdedness", "Beauty"), length = c(50, 7)) |>
    mutate(Outcome = resp) |>
    rbind(preds)

  # Data for density
  dat <- rbind(
    dat,
    data.frame(
      Typology = df$Typology,
      Crowdedness = df$Crowdedness,
      video_mean_pedcounts = df$video_mean_pedcounts,
      Predicted = df[[resp]],
      Outcome = resp
    )
  )
}
## [1] "WillingnessToWalk"
## [1] "Valence"
## [1] "Arousal"
## [1] "Restoration"
## [1] "Presence"



names(mods) <- outcome_vars
preds_psych <- preds

plots psych


dat$Outcome <- recode(dat$Outcome, "WillingnessToWalk"="Willing.-to-walk")

p_psych <- preds_psych |>
  dplyr::filter(Predicted <= 7) |>
  mutate(Beauty = fct_rev(as.factor(Beauty)),
         Outcome = recode_factor(Outcome, "WillingnessToWalk"="Willing.-to-walk"),
         Outcome = fct_relevel(Outcome, c("Willing.-to-walk", "Valence", "Arousal", "Restoration", "Interest", "Presence"))) |>
  ggplot(aes(x = Crowdedness, y = Predicted)) +
  stat_density_2d(data = dat, aes(fill = ..density..), geom = "raster", contour = FALSE) +
  geom_line(aes(size = Beauty, color = Typology, group = interaction(Beauty, Typology))) +
  scale_size_manual(values = seq(from = .05, to = 2, length.out= 7), breaks= c(1:7)) +
  scale_fill_gradientn(colours = c("white", "grey"), guide = "none") +
  scale_color_manual(values = c("urban" = col_urban, "green" = col_green)) +
  scale_alpha_ordinal(range = c(1, 0.1)) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0), breaks = c(1:7)) +
  theme_modern(axis.title.space = 10) +
  theme(
    strip.placement = "outside",
    strip.text.x = element_text(hjust = 0, size = 12),
    axis.title.y = element_text(face = "bold"),
    plot.title = element_text(face = "bold", hjust = 0.5)
  ) +
  labs(x = "Subjective Crowdedness", y = "Psychological Outcome", title = "Psychological effect of Crowdedness and its modulation by Beauty") +
  ggnewscale::new_scale_fill() +
  scale_fill_manual(values = c("urban" = col_urban, "green" = col_green)) +
  ggside::geom_xsidedensity(data=dat, aes(y=stat(density), fill = Typology), color = NA, alpha = 0.5, adjust = 2) +
  ggside::geom_ysidedensity(data=dat, aes(x=stat(density), fill = Typology), color = NA, alpha = 0.5, adjust = 2) +
  ggside::theme_ggside_void() +
  facet_wrap(~Outcome, strip.position = "left", scales = "free")

p_psych+ theme(legend.position = 'bottom') + guides(shape = guide_legend(nrow = 1, title.vjust = 1))

Table output


names(mods) <- outcome_vars


pretty_vars <- map(mods, .f = format_parameters) %>% unname() %>% unlist() 
pretty_vars
##                                         (Intercept) 
##                                       "(Intercept)" 
##                                       Typologygreen 
##                                  "Typology [green]" 
##                                              Beauty 
##                                            "Beauty" 
##                                         Familiarity 
##                                       "Familiarity" 
##                 Typologyurban:poly(Crowdedness, 2)1 
##          "Typology [urban] × poly(Crowdedness, 2)1" 
##                 Typologygreen:poly(Crowdedness, 2)1 
##          "Typology [green] × poly(Crowdedness, 2)1" 
##                 Typologyurban:poly(Crowdedness, 2)2 
##          "Typology [urban] × poly(Crowdedness, 2)2" 
##                 Typologygreen:poly(Crowdedness, 2)2 
##          "Typology [green] × poly(Crowdedness, 2)2" 
##                                Typologygreen:Beauty 
##                         "Typology [green] × Beauty" 
##          Typologyurban:poly(Crowdedness, 2)1:Beauty 
## "Typology [urban] × poly(Crowdedness, 2)1 × Beauty" 
##          Typologygreen:poly(Crowdedness, 2)1:Beauty 
## "Typology [green] × poly(Crowdedness, 2)1 × Beauty" 
##          Typologyurban:poly(Crowdedness, 2)2:Beauty 
## "Typology [urban] × poly(Crowdedness, 2)2 × Beauty" 
##          Typologygreen:poly(Crowdedness, 2)2:Beauty 
## "Typology [green] × poly(Crowdedness, 2)2 × Beauty" 
##                                         (Intercept) 
##                                       "(Intercept)" 
##                                       Typologygreen 
##                                  "Typology [green]" 
##                                              Beauty 
##                                            "Beauty" 
##                                         Familiarity 
##                                       "Familiarity" 
##                 Typologyurban:poly(Crowdedness, 2)1 
##          "Typology [urban] × poly(Crowdedness, 2)1" 
##                 Typologygreen:poly(Crowdedness, 2)1 
##          "Typology [green] × poly(Crowdedness, 2)1" 
##                 Typologyurban:poly(Crowdedness, 2)2 
##          "Typology [urban] × poly(Crowdedness, 2)2" 
##                 Typologygreen:poly(Crowdedness, 2)2 
##          "Typology [green] × poly(Crowdedness, 2)2" 
##                                Typologygreen:Beauty 
##                         "Typology [green] × Beauty" 
##          Typologyurban:poly(Crowdedness, 2)1:Beauty 
## "Typology [urban] × poly(Crowdedness, 2)1 × Beauty" 
##          Typologygreen:poly(Crowdedness, 2)1:Beauty 
## "Typology [green] × poly(Crowdedness, 2)1 × Beauty" 
##          Typologyurban:poly(Crowdedness, 2)2:Beauty 
## "Typology [urban] × poly(Crowdedness, 2)2 × Beauty" 
##          Typologygreen:poly(Crowdedness, 2)2:Beauty 
## "Typology [green] × poly(Crowdedness, 2)2 × Beauty" 
##                                         (Intercept) 
##                                       "(Intercept)" 
##                                       Typologygreen 
##                                  "Typology [green]" 
##                                              Beauty 
##                                            "Beauty" 
##                                         Familiarity 
##                                       "Familiarity" 
##                 Typologyurban:poly(Crowdedness, 2)1 
##          "Typology [urban] × poly(Crowdedness, 2)1" 
##                 Typologygreen:poly(Crowdedness, 2)1 
##          "Typology [green] × poly(Crowdedness, 2)1" 
##                 Typologyurban:poly(Crowdedness, 2)2 
##          "Typology [urban] × poly(Crowdedness, 2)2" 
##                 Typologygreen:poly(Crowdedness, 2)2 
##          "Typology [green] × poly(Crowdedness, 2)2" 
##                                Typologygreen:Beauty 
##                         "Typology [green] × Beauty" 
##          Typologyurban:poly(Crowdedness, 2)1:Beauty 
## "Typology [urban] × poly(Crowdedness, 2)1 × Beauty" 
##          Typologygreen:poly(Crowdedness, 2)1:Beauty 
## "Typology [green] × poly(Crowdedness, 2)1 × Beauty" 
##          Typologyurban:poly(Crowdedness, 2)2:Beauty 
## "Typology [urban] × poly(Crowdedness, 2)2 × Beauty" 
##          Typologygreen:poly(Crowdedness, 2)2:Beauty 
## "Typology [green] × poly(Crowdedness, 2)2 × Beauty" 
##                                         (Intercept) 
##                                       "(Intercept)" 
##                                       Typologygreen 
##                                  "Typology [green]" 
##                                              Beauty 
##                                            "Beauty" 
##                                         Familiarity 
##                                       "Familiarity" 
##                 Typologyurban:poly(Crowdedness, 2)1 
##          "Typology [urban] × poly(Crowdedness, 2)1" 
##                 Typologygreen:poly(Crowdedness, 2)1 
##          "Typology [green] × poly(Crowdedness, 2)1" 
##                 Typologyurban:poly(Crowdedness, 2)2 
##          "Typology [urban] × poly(Crowdedness, 2)2" 
##                 Typologygreen:poly(Crowdedness, 2)2 
##          "Typology [green] × poly(Crowdedness, 2)2" 
##                                Typologygreen:Beauty 
##                         "Typology [green] × Beauty" 
##          Typologyurban:poly(Crowdedness, 2)1:Beauty 
## "Typology [urban] × poly(Crowdedness, 2)1 × Beauty" 
##          Typologygreen:poly(Crowdedness, 2)1:Beauty 
## "Typology [green] × poly(Crowdedness, 2)1 × Beauty" 
##          Typologyurban:poly(Crowdedness, 2)2:Beauty 
## "Typology [urban] × poly(Crowdedness, 2)2 × Beauty" 
##          Typologygreen:poly(Crowdedness, 2)2:Beauty 
## "Typology [green] × poly(Crowdedness, 2)2 × Beauty" 
##                                         (Intercept) 
##                                       "(Intercept)" 
##                                       Typologygreen 
##                                  "Typology [green]" 
##                                              Beauty 
##                                            "Beauty" 
##                                         Familiarity 
##                                       "Familiarity" 
##                 Typologyurban:poly(Crowdedness, 2)1 
##          "Typology [urban] × poly(Crowdedness, 2)1" 
##                 Typologygreen:poly(Crowdedness, 2)1 
##          "Typology [green] × poly(Crowdedness, 2)1" 
##                 Typologyurban:poly(Crowdedness, 2)2 
##          "Typology [urban] × poly(Crowdedness, 2)2" 
##                 Typologygreen:poly(Crowdedness, 2)2 
##          "Typology [green] × poly(Crowdedness, 2)2" 
##                                Typologygreen:Beauty 
##                         "Typology [green] × Beauty" 
##          Typologyurban:poly(Crowdedness, 2)1:Beauty 
## "Typology [urban] × poly(Crowdedness, 2)1 × Beauty" 
##          Typologygreen:poly(Crowdedness, 2)1:Beauty 
## "Typology [green] × poly(Crowdedness, 2)1 × Beauty" 
##          Typologyurban:poly(Crowdedness, 2)2:Beauty 
## "Typology [urban] × poly(Crowdedness, 2)2 × Beauty" 
##          Typologygreen:poly(Crowdedness, 2)2:Beauty 
## "Typology [green] × poly(Crowdedness, 2)2 × Beauty"

md <- modelsummary(mods, coef_omit = "Cor|SD", coef_rename = rename_explanator,
             statistic = "[{conf.low}, {conf.high}]",
             estimate = "{estimate}{stars} ",
             shape = term ~ model + statistic,
             gof_map = c("nobs", "r2.conditional","r2.marginal"), fmt = 3, 
             output = "dataframe",
             ) 
md  
md <- md %>% 
  select(-part) %>%   
  mutate(term = if_else(term %in% c("Num.Obs.", "R2 Cond.", "R2 Marg."), term, pretty_vars[md$term] %>% unname())) %>% 
  mutate(term = str_replace_all(term, "video mean pedcounts", "Ped.Counts")) %>% 
  mutate(term = gsub(term, pattern = "poly(Crowdedness, 2)1", replacement = "Subj. Crowding [1st degree]", fixed = T),
         term = gsub(term, pattern =  "poly(Crowdedness, 2)2", replacement = "Subj. Crowding [2nd degree]", fixed = T)) 


md
md %>%   
   kable(format = "latex", 
        caption = 'Non-linear mixed effects models of psychological reaction.', 
        label = "models_reaction",
        booktabs = T, col.names = c("Term", rep(c("Est.", "95% CI"), length(outcome_vars)))) %>% 
  add_header_above(setNames(c(0, rep(2, length(outcome_vars))), c("", outcome_vars)) ) %>% 
  kable_classic_2(latex_options = c("striped", "scale_down")) %>%
  kableExtra::row_spec(., row = nrow(md)-3, hline_after = TRUE) %>% 
  landscape() %>%
  save_kable(file = "output/table_models_psych_reaction.tex")

joint plots


p_reg <- ((p_crowd + guides(fill = "none", color = "none") | p_beauty + theme(legend.position = "none")) / p_psych )  +plot_annotation(tag_levels = "A")

p_reg
## Warning: Removed 12 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 22 rows containing non-finite values (`stat_density()`).
## Warning: Removed 34 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 12 rows containing non-finite values (`stat_density()`).
## Warning: Removed 22 rows containing non-finite values (`stat_density()`).
## Warning: Removed 420 rows containing non-finite values
## (`new_stat_density2d()`).
## Warning: Removed 15 rows containing non-finite values (`stat_density()`).
## Warning: Removed 405 rows containing non-finite values (`stat_density()`).
## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.



ggsave(p_reg, filename = "output/plots_reg.jpg", width = 13, height = 13, dpi = 300)
## Warning: Removed 12 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 22 rows containing non-finite values (`stat_density()`).
## Warning: Removed 34 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 12 rows containing non-finite values (`stat_density()`).
## Warning: Removed 22 rows containing non-finite values (`stat_density()`).
## Warning: Removed 420 rows containing non-finite values
## (`new_stat_density2d()`).
## Warning: Removed 15 rows containing non-finite values (`stat_density()`).
## Warning: Removed 405 rows containing non-finite values (`stat_density()`).
## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.

names(dfenv) <- make.names(names(dfenv))
cor_urb <- dfenv %>% 
  filter(Typology == "Urban") %>% 
  ungroup() %>% 
  select(Willingness.to.walk, Restoration, Interest, Valence, Arousal, Crowdedness, Scenic, Width ) %>% 
  correlation(partial = TRUE, standardize_names = F) %>% 
  filter(p < 0.01) %>%
  plot()
cor_urb


cor_green <- dfenv %>% 
  filter(Typology == "Green") %>% 
  ungroup() %>% 
  select(Willingness.to.walk, Restoration, Interest, Valence, Arousal, Crowdedness, Scenic, Width ) %>% 
  correlation(partial = TRUE) %>% 
  # filter(p < 0.01) %>% 
  plot()
cor_green

Session information

sessionInfo()
## R version 4.2.2 (2022-10-31)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS 14.0
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] kableExtra_1.3.4.9000   ggraph_2.1.0            modelsummary_1.3.0.9016
##  [4] patchwork_1.1.2         see_0.7.5               report_0.5.7           
##  [7] parameters_0.21.0       performance_0.10.3      modelbased_0.8.6       
## [10] insight_0.19.1.6        effectsize_0.8.3        datawizard_0.7.1       
## [13] correlation_0.8.4       bayestestR_0.13.1       easystats_0.6.0        
## [16] ggside_0.2.2            ggdist_3.2.1            lubridate_1.9.2        
## [19] forcats_1.0.0           stringr_1.5.1           dplyr_1.1.4            
## [22] purrr_1.0.2             readr_2.1.4             tidyr_1.3.0            
## [25] tibble_3.2.1            ggplot2_3.4.4           tidyverse_2.0.0        
## 
## loaded via a namespace (and not attached):
##  [1] ggnewscale_0.4.8       minqa_1.2.6            colorspace_2.1-0      
##  [4] estimability_1.4.1     rstudioapi_0.14        farver_2.1.1          
##  [7] MatrixModels_0.5-3     graphlayouts_0.8.4     ggrepel_0.9.3         
## [10] DT_0.27                fansi_1.0.6            mvtnorm_1.1-3         
## [13] xml2_1.3.3             splines_4.2.2          cachem_1.0.8          
## [16] knitr_1.42             polyclip_1.10-4        jsonlite_1.8.8        
## [19] nloptr_2.0.3           gt_0.9.0               broom_1.0.5           
## [22] ggforce_0.4.1          compiler_4.2.2         httr_1.4.7            
## [25] emmeans_1.8.5          backports_1.4.1        Matrix_1.6-4          
## [28] fastmap_1.1.1          cli_3.6.2              tweenr_2.0.2          
## [31] htmltools_0.5.5        tools_4.2.2            lmerTest_3.1-3        
## [34] igraph_1.4.2           coda_0.19-4            gtable_0.3.4          
## [37] glue_1.7.0             tables_0.9.10          Rcpp_1.0.11           
## [40] jquerylib_0.1.4        vctrs_0.6.5            svglite_2.1.1         
## [43] nlme_3.1-162           xfun_0.39              lme4_1.1-35.1         
## [46] rvest_1.0.3            timechange_0.2.0       lifecycle_1.0.4       
## [49] MASS_7.3-59            scales_1.3.0           tidygraph_1.2.3       
## [52] BayesFactor_0.9.12-4.4 ragg_1.2.5             hms_1.1.3             
## [55] parallel_4.2.2         yaml_2.3.7             pbapply_1.7-0         
## [58] gridExtra_2.3          sass_0.4.6             stringi_1.8.3         
## [61] highr_0.10             checkmate_2.1.0        boot_1.3-28.1         
## [64] rlang_1.1.3            pkgconfig_2.0.3        systemfonts_1.0.4     
## [67] distributional_0.3.2   evaluate_0.23          lattice_0.21-8        
## [70] htmlwidgets_1.6.2      labeling_0.4.3         tidyselect_1.2.0      
## [73] magrittr_2.0.3         R6_2.5.1               generics_0.1.3        
## [76] pillar_1.9.0           withr_2.5.2            crayon_1.5.2          
## [79] utf8_1.2.4             tzdb_0.3.0             rmarkdown_2.21        
## [82] viridis_0.6.2          grid_4.2.2             data.table_1.14.10    
## [85] digest_0.6.33          webshot_0.5.4          xtable_1.8-4          
## [88] numDeriv_2016.8-1.1    textshaping_0.3.6      munsell_0.5.0         
## [91] viridisLite_0.4.2      bslib_0.4.2